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(54) Thermal imaging method and apparatus 

(57) An apparatus for dynamic imaging of the blood 
perfusion in skin measures blood flow in human or ani- 
mal skin non-invasively, by means of processing skin 
thermographic data. An area of skin is first cooled, the 
cooling means is removed, and then, using an infra-red 
detector, thermographic data conresponding to skin 
temperature is obtained at a number of points on the 



area of skin over a period of time. For each point an 
exponential fit equation for the variation with time of the 
thermographic data is calculated. The coefficients of the 
exponential fit equation are displayed in graphical nrtap 
form, and the map can be interpreted to provide infor- 
mation about the blood perfusion rate in the skin. 
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Description 

This invention relates to an apparatus and method for dynamic imaging of blood perfusion in human or animal skin, 
in particular by means of measuring and processing skin thermographic data. 
5 Many changes occur to the structure of human skin during growth. However, ignoring variation due to development 
and ageing, regional variations or the effects of pathology, there is a general pattern and structure of the vessels of the 
skin which is characteristic of most skin, for most of the time. This general pattern is presented in Fig. 1, which shows 
schematically the distribution of arteries and veins in a section of skin. 

The epidermis (epitiielial layer) 1 forms the outer waterproof boundary of the body and comprises a tiiin avascular 
10 layer less tiian 0.17 mm thick in most parts. Below this lies tiie highly vascular dermis 2, 3, 4 containing many blood 
and lymphatic vessels, where specialised structures (sweat glands, hair follicles 6, sensing organs) are located. Its 
thickness varies from 0.5 mm up to 4.0 mm. 

The dermis can be subdivided into tiie papillary layer or upper dermis 2, tiie reticular layer or mid-dermis 3 and the 
deep dermis 4. The papillary layer 2 lies close to tiie epidermis 1 and contains a vascular complex, the sub-papillary 
15 plexus, which nourishes tiie base of the epidermis. Under this Is tiie reticular dermis 3. containing a network of fibres 
and blood vessels that connect the papillary dermis 2 and the deep dermis 4. 

The deep dermis 4 contains tiie sweat glands, tiie hair follicles 6 and most importantiy, the major vascular network 
7 supplying tiie whole dermis. Below the dermis there is a layer 5 consisting of adipose tissue or subcutaneous fat, 
which is not highly vascular and whose tiiickness varies greatiy between individuals. An important characteristic of this 
20 layer is ttiat it is a poorer conductor of heat ttian other body tissues. Therefore, it acts as an insulator reducing the total 
heat lost by tiie body. 

The configuration of tiie vessels in tiie layers of the skin is shown in Rg. 2. An artery B about 100 jim in diameter 
enters the lower dermis 4. It has a lining of flat endothelial cells which lie on a relatively ttiick elastic membrane (internal 
elastic lamina). Around it there are several layers of smooth muscle, and collagen fibres. 
25 The arteries divide once or twice as they pass through the deep dermis 4 and when their branches 9 reach the mid 
dermis 3 they divide again. Here they are about 50 ^ in diameter and have only a single layer of smootii muscle cells. 
The smooth muscle layer becomes discontinuous, oblique and spiral and tiie vessels are now arterioles. No internal 
elastic lamina exists in tiie arterioles, and very few elastic fibres are detectable in tiie network supplying tiie muscle 
layer. 

30 Further on, tiie repeated subdivision of tiie arterioles reduces tiieir diameter down to approximately 10-15 ^m form- 
ing the pattern of the capillaries 10 in the upper dermis (papillary dermis) 2. This is a network of horizontal capillaries 
lying parallel to the surface (sub-papillary plexus). Above it are tiie terminal capillaries 11. These are vertical loops 
which arise from the arterioles and drain it into tiie horizontal sub-papillary venus plexus. 

The venules 1 2 in tiie upper and middle dermis are more numerous than the arterioles. Their diameter ranges from 
35 40-60 ^m in the upper and mid-dermis to 100-400 ^m in the deeper tissues. 

Some 80% or more of total body heat transfer occurs through the skin. The remaining part takes place in mucous 
memt>ranes. Heat can be lost from the skin by four different means, namely radiation, conduction, con vention and 
evaporation. 

The basic mechanism of heat exchange in tiie human body through skin to the outside environment 24 can be seen 
40 in Fig. 3. When homeostasis requires that the body 20 conserves heat, blood flow in ttie skin decreases by vasocon- 
striction 21 to limit heat loss. When heat must be lost from the body 22 to maintain the stability of the internal environ- 
ment, flow of warm blood to the skin increases by vasodilation 23. 

The temperature of the skin surface is lower than that of the body core 25 because there is a constant heat loss 
from tiie body surface to the environment. 
45 Conduction from the underlying tissues provides continuous, relatively small input of heat to the skin, but the main 
mechanism of heat exchange to the skin is by perfusion of blood, at core body temperature, from the main circulation. 
Additionally, for reasons associated witii tiie basic thermal physics, heat transfer from the warm arterial blood to tiie sur- 
rounding tissues 26 occurs mainly in vessels greater than 50 ^m in diameter. 

For tiiis reason skin temperature, especially at high perfusion rates, reflects heat exchange at tiie deeper levels of 
50 the dermis. Consequentiy an area of tiie skin witii relatively high blood perfusion rate has a surface temperature higher 
than tiiat of an area with lower perfusion. 

Many areas of clinical medicine would benefit from a metiiod of quantifying regional perfusion to tiie skin, for exam- 
ple, plastic surgery, amputation surgery, and cancer diagnosis. Different techniques have been developed tiirough the 
years for this purpose, but most skin blood flow measurement techniques suffer from a number of limitations. 
55 Radioisotope ciearance involves measurement by a radiation detector of tiie rate of removal of an intradermally 
injected or transcutaneously absorbed amount of radioisotope. Flow can be quantified by this method but measure- 
ments are affected by needle trauma or by fat solubility problems in tiie case of Xenon^^. 

Optical metiiods such as Laser Doppler Fiowmetry and Photoplethysmograptiy, altiiough very popular because of 
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their minimal interference with the skin, have some considerable disadvantages too. 

Laser Doppler Flowmetry uses the Doppler shift of light reflected from a small volume of illuminated skin tissue to 
indicate the flux of moving red cells and therefore dermal perfusion in that portion of the skin. It has, though, uncertain 
penetration depth, measures in very snnall volumes of tissue, has a poor response to movement artefact, and gives only 
5 relative values of flow. Instrumentation applying this technique has been developed, using a scanning mirror rather than 
fibre optics in order to guide the laser beam. Using a scanning min^or may solve the spatial heterogeneity problem but 
it is still under trial and the produced measurements are only qualitative and slow to obtain. 

Photoplethysmography uses visible or infra-red radiation reflected from the skin tissue to assess the blood flow in 
the skin layers. Although this technique indicates the timing of events, it provides a poor measure of change in volume. 
10 and it is very sensitive to motion artefact. 

The thermal clearance method relates artificially induced temperature gradients longitudinal to a small region of the 
skin to blood flow in that region. This is also a point measurement technique which involves changes in dermal per- 
fusion. 

Medical thermography is a diagnostic method of imaging the distribution of temperature of human skin, tt has a 
IS number of distinct advantages, tt is completely non-invasive, a real time procedure, and there is no radiation hazard. 
This method can form the basis for regional mapping of skin blood flow, relating skin temperature with perfusion. How- 
ever, static thermographic imaging is very dependent on environmental changes of temperature, the thermoregulatory 
status of the subject, and mechanisms of heat exchange that are not related to blood perfusion. 

It is an object of the present invention to provide an apparatus and method for dynamic imaging of the blood per- 
20 fusion rate in skin, which overcomes the problems of the prior art described above. 

According to a first aspect of the present invention there is provided a method of measuring blood perfusion in 
human or animal skin, comprising the steps of: 

(a) applying a cooling means to an area of skin; 
25 (b) removing the cooling means from the area of skin; 

(c) using an infra-red detector to obtain a first set of thermographic data TD^ corresponding to skin tenrperature at 
a plurality of points on the area of skin at time t^ ; 

(d) repeating step (c) at least once to obtain a total of n sets of tiiermographic data TD^ to TD^ corresponding to 
times ti to tp, where n is at least 2; 

30 (e) calculating tiie coefficients of an exponential fit equation for the variation with time of tiie thermographic data 
TDi to TDn at each of the plurality of points; and 

(f) displaying the coefficients of the exponential fit equation in graphical form. 

Preferably the coefficients of the exponential fit equation are calculated using a curve stripping technique. Prefera- 
35 biy n is at least 4, most preferably at least 10. 

Preferably the coefficients of the exponential fit equation are displayed using colour mapping. 
Preferably the cooling means is a cold object which s placed on tiie skin. 

Preferably tiie method furtiier includes tiie step of placing a reference object on the area of skin so that tiie thermo- 
graphic data TDp is obtained at the same points on the area of skin at each successive repeat of step (c). The reference 
40 object may be a reflective strip. 

According to a second aspect of the present invention there is provided a method of processing and displaying data 
relating to blood perfusion in human or animal skin, comprising the steps of: 

(a) generating a series of n thermographic image matrices IM^ to IM^ corresponding to times ti to tn, each thernno- 
45 graphic image matrix comprising temperature data relating to a plurality of points on an area of skin; 

(b) creating a filter matrix corresponding to a region of interest selected from ttie area of skin; 

(c) applying tiie filter matrix to each tiiermographic image matrix IM^ to IMp to obtain an cropped thermographic 
image matrix IMCRi to IMCRn having the value zero for points outside tiie region of interest; 

(d) constructing an array AR^ for each point m in tiie region of interest, the array containing temperature data T^ 
50 at ttiat point corresponding to times t^ to tp; 

(e) for each an^ay ARjj^ performing a regression analysis to calculate the coefficients of the exponential components 
of tiie best fit equation 



55 



T, = T„„ + Ae'" + Be"*Ce'=' 
where A, a, B. b. C. c are coefficients and Tf^jn is the mininnum temperature data in each array; 
(f) constructing a coefficient matrix (Z^, Zq etc.) having as elements a parameter derived from the coefficients (A. 
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a etc.) calculated In step (e) for each point m in the region of interest; and 
(g) displaying the coefficient nnatrix as an image. 

Preferably the parameter in step (f) is equal to one of the coefficients (A, a etc.) calculated in step (e). Alternatively 
5 the parameter in step (f) may be equal to the sum of the coefficients A, B and C calculated in step (e), representing the 
total tissue volume. Alternatively the parameter in step (f) may be equal to tiie sum of tiie coefficient products (aA + bB 
+ cC) divided by tiie sum of the coefficients (A + B + C), representing the total flow. 

Preferably the region of interest corresponds to an area of skin which has previously been cooled. 

Preferably the coefficients of the exponential components are calculated using a curve stripping technique. 
10 Preferably the minimum temperature data Imm ^<^^ ai^ay 'S subtracted from each temperature data in the array 
before calculation of the coefficients. 

Preferably tiie data corresponding to a first interval of time after cooling Is discarded before calculating coefficients. 
The interval of time is preferably less than 10 seconds, more preferably less tiian 5 seconds, most preferably 3 or 4 sec- 
onds. 

15 According to a third aspect of the present invention there is provided an apparatus for measuring blood perfusion 
in human or animal skin, comprising: 

(a) an infra-red detecting means adapted to obtain periodically a set of tiiermographic data (TD^ to TDn) corre- 
sponding to skin temperature at a plurality of points on an area of skin; 
20 (b) data storage means adapted to store each set of thermographic data (TD^ to TDp); 

(c) computational means for calculating the coefficients of an exponential fit equation for tiie variation with time of 
the thermographic data (TDi to TDJ at each of the plurality of points; and 

(d) image display means adapted to display the coefficients of the exponential fit equation in graphical form. 

25 The method of the present inverrtion makes use of tiiermographic measurements (images) of dynamic changes in 
skin temperature, i.e. changes tiiat are a result of the reheating of the skin after contact between the skin and a cold 
object. 

The acquired images may be processed by suitable software. The purpose of this process is to perform a multi- 
exponential regression analysis to data that form reheat curves. The resultant coefficients of the components may then 
30 be reconstructed as new images. 

The imaging techniques of tiie present invention may be applied to breast tumours, critical limb ischaemia and dia- 
betic ulcers. 

An embodiment of tiie invention will now be described with reference to the accompanying figures where: 

35 Fig. 1 is a schematic drawing of tiie distribution of arteries and veins in the skin; 
Fig. 2 shows the most common vascular pattern and structure in adutt skin; 
Fig. 3 is a simplified diagram showing the skin as a thermo-regulatory organ; 
Fig. 4 shows black body spectral radiant omittance (Planck's Law); 

Fig. 5 illustrates Planck's Law as a semi-log plot, showing a curve family from 100K to 1000 K; 
40 Fig. 6 illustrates tiie operating principles of a pyroelecti-ic detector; 
Fig. 7 is a diagram to illustrate tiie concept of hot "core"; 

Rgs. 8 to 10 show examples of Curve Stripping (step I, A component, step II. B component and step III, C compo- 
nent respectively); 

Fig. 1 1 is a typical screen image showing the reheat data curve at a cfiosen location and the coefficients calculated 
45 by the software curve stripping method for the reheat data curve; 

Fig. 12 shows a confrol tiiermographic image of an area on a forearm for normal skin data; 

Rgs. 13 to 15 shows the images produced by the method of tiie invention for ttie coeffiderrts A and a; B and b; and 

C and c respectively for normal skin data; 

Fig. 16 shows a control tiiermographic image of an area on a forearm for tuberculin reaction skin data; and 
50 Figs. 1 7 to 1 9 shows the images produced by the metiiod of tiie invention for tiie coefficients A and a; B and b; and 
C and c respectively for tuberculin reaction skin data. 

The invention makes use of infra-red thermography, a known technique used to monitor the temperature disti'ibu- 
tion of a material by detecting the radiation that this material emits. The electi^omagnetic spectrum is divided Into a 
55 number of wavelength regions, called "bands", distinguished by the metiiods used to produce or detect the radiation. 
The only difference between the bands is in wavelength. All tiie physical laws that govern these different regions are 
exactly the same. The background theory to the method of the Invention is described below. 

The infra-red band lies between 0.75 ^m (the end of visual band) to 100 ^m (tiie t>eginning of the microwave band). 
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It is subdivided into four sections, namely near infra-red (0.75 - 3 fim), middle Infra-red (3-6 }im), far Infra-red (6-15 )im) 
and extreme infra-red (15-100 ftm). 

Electromagnetic radiation is emitted spontaneously by all objects having a temperature above absolute zero (- 
273*'C or OK). 

A blackbody is defined as an object which absorbs all incident radiation, and emits all its thermal energy in the form 
of electromagnetic radiation, at any wavelength. In other words a blackbody is a theoretical concept, that behaves as a 
perfect radiation absorber and emitter. 

Human skin, at Zl'^C, has an emissivity of 0.98, and for this reason is a structure very dose to a blackbody. 
The spectral distribution of the radiation from a blackbody is described by means of the following formula: 

W^b = -™^x10-«[WattsMi%m] 

X (e -1) 

15 where: 

WAt, = the blackbody spectral radiant emittance at wavelength X 
c = the velocity of light (3 x 1 0 ® m/sec) 
h = Planck's constant (6.6 x 10'^^ Joule sec) 
20 k = Boltzmann's constant (1 .4 x 1 0'^^ Joule/K) 
T = the absolute temperature of the blackbody (K) 
X = wavelength (m) 

This formula describes the relationship between the radiation energy and wavelength that is emitted by a blackbody 
25 in a certain temperature. Plank's formula, when plotted graphically for various terrperatures, produces a family of 
curves. Following any particular Planck curve the spectral emittance is zero at X = 0. then increases rapidly to a maxi- 
mum at a wavelength X^a^, and then approaches zero again at very long wavelengths. The higher the temperature, the 
shorter the wavelength at which the maximum occurs. Rg. 4 shows one family of ounces for a typical body. Spectral radi- 
ant emittance (in watts/cnrt^) is plotted against wavelength (in micrometres) for a series of temperatures. 
30 By differentiating Planck's formula with respect to X, and finding the maximum. Wien's formula is derived: 

2898 r 1 
>-max = — Q^m] 

35 

This formula expresses mathematically the common observation that colours vary as the temperature of a thermal 
radiator increases. For a terrtperature close to 32°C (about 305K). which is close to the human skin tenrperature, X^nax 
= 10 fim. This is the wavelength of the radiation of highest intensity that the skin emits, and it belongs to the far infra- 
red region. 

40 In Fig. 5 a semi-log plot of the Planckian curves is presented, and the part of the spectrum that corresponds with 
the radiation emitted by human skin appears (300 K). Spectral radiant emittance (in watts/cm^ is plotted against wave- 
length (in miCTometres) for a series of temperatures. 

By integrating Planck's formula from X = 0 to X = ~, the total radiant emittance (Wb) of a blackbody is obtained. 

45 Wb = ar^[Watts/m^] 

where a = Stefan-Boltzmann constant (5.7 x 10"* W/m^ 

This is the Stefan-Boltzmann formula, which states that the total emissive power of a blacklxxiy is proportional to 
50 the fourth power of its absolute temperature. Graphically. Wb represents the area under the Planck curve for a particular 
temperature (see Figs. 4 and 5). 

Infra-red (IR) detectors are available in many different forms, but only some of these meet the requirements for 
medical thermal imaging. These have a fast response time, are most sensitive at X between 1-20 \im, and (for some 
devices) have an ability to be incorporated into anrays of detectors. IB detectors are used to convert infrared energy into 
55 electrical signals, and can be either thermal detectors or photon detectors. Devices which detect infra-red radiation 
remotely range from simple point bolometers and thermoelectric detectors to the more sophisticated imaging systems 
such as scanners and pyroelectric vidicons. The method of the invention may make use of any of these devices, which 
are described in more detail below. 
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Pyroelectric detectors (see Fig. 6) are made from ferroelectric crystals which have an electrical polarisation 
strongly dependent on temperature. 

A temperature difference (Delta T) Is induced between the thermal insulation layers 31 and the absorbing black 
layer 32. This is proportional to the chopped incident radiation power 33 (and wavelength] absorbed by the black layer 
5 . A voltage signal V is then generated by the ferroelectric crystal 34 of pyroelectric rriaterial between the electrodes 35, 
36, and this signal is also proportional to the temperature difference and therefore to the inciderrt radiation wavelength. 

Bolometers comprise thin metallic resistance elements onto which radiation is focused by an infra-red lens or mir- 
ror. Their resistance varies with temperature. Both devices operate using tiiermal detectors. 

An electrical output is obtained from these devices by mechanically chopping the radiation falling on the crystal. 
10 Botii devices are portable and can be used to give accurate temperature measurements by pointing them at a small 
area of skin, usually from close range. 

Photoconductive detectors such as Indium Antimonide (InSb) or Mercury Cadmium Telluride (MOT) are photon 
detectors which operate on tiie principle ttiat the material's crystal conductivity varies according to the number of inci- 
dent infrared photons in its sensitive waveband. These detectors require cooling by liquid nitrogen which limits tiieir 
15 mobility but by using a number of focusing and scanning mirrors they have considerable flexibility. 

The Thermographic Imaging system used most successfully in the method of the invention is of the type described 
above (Agema Thermovision 880). The detector used is an MOT element (long wave version) covering 8-12 fim of tiie 
infra-red bank (Agema operating manual). In such a scanner electromagnetic energy radiating from tiie object being 
scanned is focused by an infra-red lens into a mirror. The mirror is oscillated by a dc motor. The optical output from the 
20 oscillating mirror is focused by three fixed mirrors onto a horizontal mirror polygon which rotates at 1 5000 rpm. Botii ttie 
oscillator mirror and tiie horizontal mirror polygon are controlled by a microprocessor. The miaoprocessor provides hor- 
izontal and vertical trigger pulses to the Control Unit trigger circuit. The temperature resolution of the scanner is 0.05'*C. 

This scanner operates across the temperature range 15-40**C and tiie incoming signal from the detector is digi- 
tised. After that the tiiermographic image is displayed on a colour monitor, and the digitised data can be stored on com- 
25 puter discs in order to be processed later. 

The Pyroelectric Vidicon tube is a newer development in tiiermographic imaging and is also suitable for use in the 
method of the invention. It operates in a similar manner to tiie optical vidicon camera with the exception that the target 
material in the tube is pyroelectric (ferroelectric crystals whose electrical polarisation alters with temperature). 

There is an important requirement for all the infra-red imaging systems capable of precise temperature measure- 
so ment. A tiiermal reference is required for calibration. This can either be used in tiie field of view (near to tiie material to 
be scanned) or Incorporated into the optical system. 

The method of the invention makes use of dynamic tiiermography, which is concerned with the thermographic 
imaging of tiie human body by detecting the infra-red radiation tiiat is emitted by tiie skin. In an ambient temperature of 
22°C a healthy person might have a skin temperature that ranges from 35°C over the sternum to 25''C over the feet. As 
35 discussed before, at these temperatures peak infra-red emission occurs around 10 ^m wavelengths. Using a scanner 
that is sensitive in this region of radiation, a temperature map of tiie skin is acquired. Interpreting tiiis information, an 
experienced physician can diagnose different sorts of abnormalities. A simple but useful concept when discussing inter- 
nal body temperature is that of the cenfre "core" (see Fig. 7). 

This "core" 40 includes the vital organs of the body 41 located in tiie trunk and head. Its mean temperature, for 
40 healthy subjects at rest, is around 37**C. The study of mechanisms by which the thermal energy is transferred from 
deeper tissues to tiie skin for subsequent exchange with the environment can indicate much about physiological mech- 
anisms in health and disease. 

Infrared thermography can be used to detect changes in the vascularity of the skin caused by a variety of Internal 
disorders. Static thermography has been successfully applied in tiie tiiermographic assessments of joints, particularly 
45 for obsen/ing the response to ti-eatment of conditions such as synovitis. Static thermography has been used for breast 
tumour diagnosis but failed in favour of mammography, since it pointed to many false positive diagnoses. Today static 
thermography is used for limb amputation level assessments, diagnoses of deep vain tiironrtooses, and inflammations 
(e.g. back pains). However the simple and quick production of an Image which indicates how terrperature changes 
across an area of skin and which can be readily interpreted by medical practitioners or surgeons would be of great ben- 
50 ef it. and would provide a further non-invasive diagnostic tool. 

The following is an example of a dynamic tiiermographic procedure according to tiie invention. 

A cooling element (25''C) was placed on the subjects forearm for 15s so that the skin temperature decreased. 
Immediately after removal of the cooling block, tiie Image capture program was activated v^en tiie skin site was 
55 held in a stable position. 

Thirty images were collected: 10 images at Is intervals, 4 images at 3s. 8 images at 1 0s, 5 images at 20s, 3 images 
at 30 sec. The total time of acquisition was five minutes. 

The collected data were converted to ASCII format using a conversion program. During the conversion tiie option 
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that produced images having the higher pixel resolution was selected (136 by 136). 

The above procedure was carried out in two series of images, one recorded from normal skin and the other from 
inflamed skin. 48 hours after intradermal injection of tuberculin Purified Protein Derivative (PPD). 

5 Coefficients of an exponential fit equation for the variation with time of the thermographic data were calculated 
using software deconrposftion techniques. 

Fitting curves to multi-exponentiat data is a challenging problem encountered in many medical research areas 
fiekis such as magnetic resonance imaging (MRI) studies, pharmacokinetics, radioactive tracer techniques, blood flow 
studies and H2 clearance studies. The solution of the present invention is applicable to these problems also. 

10 According to the invention a multi-exponential curve deconposition technkiue was applied. In order to fit a curve to 
measured data. This data had been acquired during dynamic thermographic measurements of human skin as 
desaibed above. The data was treated with a regression method called curve stripping, curve peeling or successive 
subtraction. In a graphic solution data was plotted on semi-log graph paper. A line was then drawn through the last 
points of the curve ("tail"), and its extrapolated values were subtracted from the first points of the curve ("fronfO. The 

IS process was repeated until there was no n^ore data. 

The approach can be Illustrated graphically. Fig. 8 illustrates the first step in the curve stripping technique and 
shows one set of log-data. Measured data is plotted logarithmically on the Y axis against the varied parameter on the 
X axis. A line 50 is drawn from the last point 51 towards the first points. This line is the first fit-line and has a satisfactory 
con^elation with the last 8 points 51 through to 52. This line gives the A component of the exponential equation. 

20 Fig. 9 illustrates the next step, the calculation of the B component of the exponential equation. After the calculation 
of the fitted values for the first tine, the method continues subtracting these values from the original (experimental) ones. 
The residuals are the data set for the second fit-line, again plotted on a logarithmic Y axis against a linear X axis for the 
varied parameter. Another line 53 is drawn from the end 54 towards the front and this time has a satisfactory correlation 
with 1 1 data points 54 through to 55. This line gives the B component. 

25 It should be noted that in the graphical method the selection of data points which either belong or do not belong to 
a fit-line is done by visual observation of the data set. The first 7 points that are left for the next fit-line seem to belong 
to a different data "family" from those data points which fit with this second fit-line Although one line could be fitted to 
all the points of Fig. 9 it seems that this would not be right. This Visual" nature of the curve stripping method is dis- 
cussed in more detail below. 

30 The last step of this fitting procedure is shown In Fig. 10, which shows the calculation of the C component of the 
exponential equation using fit-line 56 In a similar way to that of the B component described above. 

After the calculation of the 3 fit-tines (A, B and C components), their exponentials are derived and the multi-expo- 
nential fit equation is the summation of these exponentials. 

Although the curve stripping technique is simple it gives satisfactory results for many different types of experimental 
35 data. It requires a minimal amount of computer processing capacity, so that there is no reason to use a method that 
requires access to a large computer, a lot of computing time and which produces, in the end, only a slightly better fit. 

An advantage of the curve stripping technique is that it is a visual technique, requiring the user to trust one's eyes 
rather than numbers. Numerical calculation may be more accurate, but the experienced eye can "detect", "evaluate" 
and finally "decide" much more eff idently and correctly than any "sophisticated" mathematical procedure. 
40 The visual quality of this method has been of major interest in the implementation of the software and much con- 
sideration was given into how best to include this feature in the software. 

Fig. 1 1 shows a typical screen image obtained from the software showing the reheat data cun^e at a chosen loca- 
tion and the coefficients calculated by the software curve stripping method for the reheat data curve. The screen image 
shows the cun/es used to calculate the coefficients as well as the coefficients themselves. The three upper lines show 
45 the value (square of the con-elation coefficient) for coefficients B. C and A respectively from the top. The reheat data 
curve is indicated by tiie lower line from 26''C to 35''C. The log-data curves are shown by the intermediate four lines. 

In order to create the algorithms for reconstructing the parameters, various software was developed. This will now 
be described In more detail. The software programming languages QBASIC, Turbo C++ and TurboBASIC were used. 
Two sets of data were processed using the software, one from normal skin thermographic images, and one from 
50 Images of the tuberculin reaction. Six basic procedures were implemented by the software as follows: 

Conversion from AGEMA image file format to either a binary or ASCII file format. 
Creation of a Large Image File. 
Image cropping using a polygon mask. 
55 • Creation of a Reheat Data File. 

Curve stripping or curve decomposition into exponential components. 
Parametric file image display. 
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Conversion from AGEMA im age file format to either a binary or ASCII file format 



The AGEMA thermographic camera produces its own particular image file format containing the tenperature infor- 
mation (each pixel of the 136*136 matrix represents a temperature value) and calibration data. Programs CONVERT1 

5 and C0NVERT2 were used to convert the original AGEMA inrtage files into ASCII and binary files respectively. The 
C0NVERT2 program also produced an additional file that contains tiie paths and file names of the original AGEM Afiles 
and tiie converted output binary temperature files. This filename is automatically called by the output file name 
"FILENAME" which has been generated by the user and has a .VAR extension and is located in tiie same directory as 
the one specified by the user when entering tiie output filenames and directory details. The file also contains the timing 

10 information which is required for the curve stripping algorithm. The format for tiie FILENAME.VAR file is as follows. 

NUMBER„OF_FILES 

MINTEMP 

MAXTEMP 

15 "SOURCE_PATH_AND_FILENAME.001 "/TARGET_PATH_AND_nLENAME.001 " 
"SOURCE_PATH_AND_RLENAME.003",'TARGET_PATH_AND_FILENAME.002" 

Etc... 

The NUMBER_OF_FILES number specifies the number of files that have been converted by the AGEMA conver- 
ge sion utility and tiie MINTEMP and MAXTEMP numbers specify the minimum and maximum temperature range for tiie 
output binary temperature data files specified by tiie TARGET_PATH_AND_FILENAME string. The timing information 
is held in the extension of the SOURCE_PATH_AND_FILENAME character string. 

The binary temperature data has a resolution of 254 discrete tenperatijre values from tiie minimum temperature 
(normally 24 degrees Celsius binary value 00000000 to 37 degrees Celsius binary value 11111110). The value 255 
25 binary 11111111 is not used as tiiis is required for the image cropping algorithm. 

Creation of a Large Image File 

The software developed for the analysis of the dynamic tiiermographic image processing is called PR0JECT1. 
30 The program PROJECT1 was developed as an easy to use program with a simple menu system allowing the user 
to analyse and display tiie thermographic images and to quickly produce tiie parametric image files. After the conver- 
sion of tiie original AGEMA files to the output binary temperature image files using tiie AGEMA conversion program 
C0NVERT2 (which may be run from a main menu available to tiie user), a composite file called a Large Image File is 
created, called FILENAME.LIR 

35 The FILENAME.LIR file Is a binary file which combines the binary image files sequentially to produce one large file 
which is easier to manipulate using a single file pointer Instead of opening dozens of individual files. This improves the 
speed and efficiency of the dynamic tiiermographic image processing software and saves the program PROJECT1 
opening dozens of individual binary Image files. The creation of the Large Image File is achieved by opening a file called 
FILENAME.LIF in the same directory as tiie binary image files. Each binary image data file is tiien opened sequentially 

40 and tiie data is added to the output FILENAME.LIF file stream until there are no more binary image files. Therefore, if 
there are 30 binary image files consisting of 1 36x1 36 pixel matrices ( 1 8496 bytes) the FILENAME. LIF file that would be 
generated is 30x18496 bytes in length, i.e. 554881 bytes in length. (Note the extra byte marks the End of File or EOF). 
A single file now contains all the thermographic image information. The next step in tiie parameti-ic image processing 
is to define a Region of Interest (ROI) and to generate and store the reheat temperature data sets for each pixel in a 

45 separate file. The program PR0JECT1 also has an algorithm to deal with tiie problem of slight but noticeable move- 
ments of the experimental subject during the temperature data acquisition. In order to detect any movement the user 
must specify an option (F8-PIXEL TEMPERATURE DISPLAY/APPLY MOVEMENT DETECTION) before the creation of 
the Large Image File. 

50 Image Cropping Using a Polygon Mask 

The dynamic thermographic image processing metiiod obtains more accurate infornnation about blood perfusion in 
the skin tiian can be achieved using only a single static tiiernwgraphic image. The only region of skin which is in a 
dynamic mode is the area of skin which has been subjected to tiie cold challenge (i.e. has had a cooling element 
55 applied thereto). The area of tiie cold challenge is easily observed in tiie first tiiermographic image by the distinctive 
rectangular cooled area of skin where the cooling element has been placed prior to the acquisition of tiie series of the 
thermographic images. The remaining static regions of the thermographic images are ttierefore excluded to provide a 
region of interest (ROf) by applying a polygon mask and setting tiie excluded temperature data values to 255. In prac- 
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tice this was achieved by using the QBASIC command POINT When the user is prompted to create the data reheat file 
the program PR0JECT1 scans the thermographic image polygon mask using the POINT command and discards any 
reheat temperature data sets where the value of the pixel which Is being scanned is set to 255. Therefore pixels within 
the region specified by the polygon described by the user are allowed to be used in the creation of the reheat data file 

5 (FILENAME. BHD). The reheat data file is in an ASCII format which contains the timing data. Information text header, 
and the pixel spatial location X,Y followed by the reheat temperature data set. The information text header consists of 
a number of text strings which allow the user to enter information about the experimental subject such as name, age. 
location of the cold challenge and address etc. This is done by selecting a first option (F12-CHANGE INTERNAL SET- 
TINGS/EDIT TIMING SEQUENCE), then selecting a sub-option (4-EDIT SUBJECT DETAILS) and then entering the 

JO appropriate text at each prompt. This option should be carried out before the creation of a Reheat Data File or the text 
header information will be missing. 

QrgatiQn Qf a Rgh^^t Data File 

IS The reheat data file is generated by the program PR0JECT2 and contains the following information in the following 
ASCII format. 

"PATH_AND_FILENAME_OF_PARAMETRIC_DATA_RLE" 

"DATE_OF_RLENAME,RHD_CREATION" 
20 "nME_OF_RLENAME.RHD_CREATION" 

"HEADER_TEXT_FIRST_LINE" 

"HEADER_TEXT_SECOND_UNE" 

"HEADER_TEXT_THIRD_LINE" 

"HEADER_TEXT_FOURTH_LINE" 
25 "HEADER_TEXT_RFTH_LINE" 

"HEADER_TEXT_SIXTH_UNE" 

"HEADER_TEXT_SEVENTH_LINE" 

NUMBER_OF_RLES 

TIME_VALUE(1) 
30 TIME_VALUE(2) 

TIME_VALUE(3) 



TIME_VALUE(NUMBER_OF_RLE) 

XL0CATI0N,YL0CATI0NJEMPERATURE(1) 

35 TEMPERATURE(NUMBER_OF_RLES) etc 

Curve Stripping or Curve Decomposition Into Exponential Components 

The next and most important step is to implement tiie algorithm to perform the multi-exponential curve deconnpo- 
40 sition using a modified regression analysis technique. The software is used to perform a multi-exponential curve fit in 
order to define three exponential component functions for each reheat temperature data set. Each reheat temperature 
data set corresponds to one individual pixel within tiie specified Region of Interest (ROI) defined by the user. The object 
of the software is to construct a number of new images using the parameters derived from the exponential decomposi- 
tion (hence the term parametric imaging). The parameters are tiie coefficients A,a,BiC,C,c of the mathematical model: 

45 

Tr (t) sminval + Aexp'*^+Bexp"***'+Cexp'*'' 



50 

where lt(X) is tiie reheat temperature variable, 

where minval is the first temperature value in the reheat temperature data set, and 
where A,afB,c,C,c are the coefficients of the 3 exponential components. 

55 The process of generating the coefficients A,a,B,c,C,c is a difficult one and the following metiiod assigns values to 
the six variables in order to generate a reasonable con'elated lit" between the experimental reheat temperature data 
set and tiie values which are obtained when tiie six coefficients are used to reconstruct tiie theoretical or "predicted" 
temperature values using the above mathematical model. 
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The program used to perform the multi-exponential curve deconposition Is called CURVES. 
The Reheat Data File (FILENAME. RHD) contains all the information that is required for the program CURVES to 
operate. The program CURVES performs the following tasks: 

The program reads the first 10 lines of the file which contain the following: 

"PATH_AND_FILENAME_OP_PARAMETRIC_DATA_PILE" ' dataf ile$ 
"DATE OF FILENAME. RHD CREATION" 'datel$ 



"TIMEOF^FILENAME . RHDCREATION" 


' timel$ 


"HKADERTEXTFIRST^LINE" 


' datatext$ 


"HEADERTEXTSECONDLINE" 


' datatextname$ 


"HEADERTEXTTHIRDLINE" 


' datatextage$ 


"HEADER_TEXT_FOURTH_LINE» 


' datatextlocation$ 


"HEADERTEXTFIFTHLINE " 


'datatextaddressl$ 


"HEADER_TEXT_SIXTH_LINE " 


' datatextaddress2$ 


" HEADERTEXTSEVENTHIilNE " 


' datatextaddressB $ 



The program then reads the following data: 

NDMBEROFPILES 
TIME_VALXJE(1) 
TIME_VALUE(2) 
TIME VALUE (3) 



TIME VALUE (NUMBER OF FILES) 



' size% 
'timeval%(l) 
' tiineval% (2) 
' timeval% (3) 



' timval% (size%) 



The NUMBER_OF_FILES determines how many TIME-VALUE data values are to be read by the program. The 
TIME_VALUE data values are then read by the program and are then stored into the array timeval% (1-size%) array. 
The working arrays required for the curve decomposition are then dimensioned using the NUMBER_OF_FILES integer. 
This variable is called size% in the implemented program. 

The output parametric data file FILENAME.DAT is opened in the directory specified in the 
PATH_AND_FILENAME_OF_PARAMETRIC_DATA_FILE (datafileS) string and the following header text output is writ- 
ten to the file: 



"DATE OF FILENAME . RHD_CREATION" ' datel$ 

"TIME OF FILENAME . RHD_CREATION" ' t imel$ 

"HEADER_TEXT_FIRST_LINE» ' datatext$ 

"HEADERTEXTSECONDLINE" ' datatextname$ 

"HEADERTEXTTHIRDLINE" ' datatextage$ 
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"HEADER_TEXT_FOURTH_LINE" ' datatext locat ion$ 

»HEADER_TEXT_FIFTH_LINE" ' datatextaddressl$ 

''HEADER_TEXT_SIXTH_LINE" ' datatextaddress2$ 

"HEADER TEXT SEVENTH LINE" ' datatextaddress3$ 



10 

The program then performs the following loop: 

The program reads the XLOCATION (xx%), YLOCAT10N (yy%) and then the reheat temperature data set TEM- 
PERATURE(1)...TEMPERATURE(NUMBER_0F_FILES) In to the array trhv (1..size%). 

The program CURVES then carries out the subroutine called curvestrip. The subroutine curvestrip calculates the 
IS coefficients A,a,B,c,C,c and these parameters together with the XLOCATION and YLOCATION are then printed to the 
screen and written to the output data file FILENAME.DAT using the ASCII format. 
XL0CAT10N,YL0CATI0N,A,a,B,c,C,C 

The program then terminates when there are no more reheat temperature data sets, Le, the end of the 
FILENAME.RHD file. The FILENAME.RHD and FILENAME.DAT files are then closed. 
20 The FILENAME.DAT file is referred to as a parametric image file and the data contained In the file is used by the 
program RESULTS1 which will image each parameter to determine the nature of the spatial variations of each paramet- 
ric image. 

The subroutine curvestrip performs the following algorithm to produce the required coefficient parameters 
A,a,B,c,C,c. 

25 The method of cun/e deconnposition used in the software is very sensitive to any offsets which occur in the analysis 
of the reheat temperature data set. Therefore the minimum value is subtracted from each temperature value. The vari- 
able offsetemp is set to the first temperature data point trhv(1). Assuming that the proceeding temperature data points 
are incremental in value the variable offsetenrp is subtracted from the reheat temperature data set held in the array trhv 
(2..size%) to produce the temporary variable diff. If the variable diff is negative, i.e. a temperature data point has a tem- 

30 perature which is less than the minimum temperature tiien the diff variable is set to zero. Since the natural logarithm of 
zero is undefined the diff variable is incremented by adding 1 . Therefore any temperature less than the first tenperature 
data point has the value of the natural logarithm of 1 , i.e. zero. The natural logarithnre of the reheat temperature data 
set are stored in the an^ay lgrhv(1 ..size%) with tiie value 1 lgrhv(1) always set to zero. 

The curve decomposition begins its operation by working from the last data point and works backwards towards the 

35 first data point. In order to simplrfy the procedure it was decided to reverse the order of the two arrays 
timevall%(1..size%) and lgrhv(1..size%) so that the algorithm can begin its operation by starting at the beginning of the 
arrays timeval%(1 ..size%) and lgrhv(1 ..size%) and therefore increment the array pointers in a positive manner. This is 
achieved using the subroutine reverseorder. The an-ays called temporary(1 ..size%) and temporary%(1 .size%) are used 
as temporary storage arrays to allow the order reversal in the an'ays timeval(1 ..size%) and lgrhv(1..size%). 

40 The evaluation of the first line can now take place. 

The data that are to be processed are the natural logarithms of the original data. This means that a linear regres- 
sion equation will regress the original data exponentially. 



45 



50 



Assume the regression equation: Y«at+b 

Calculating tiie exponential of this equation, the following equation is derived: 

exp^ = exp*"****^ = exp*^ x exp** = A x exp*^ 



55 In order to evaluate which is tiie best fit for tiie logarithmic data the conrelation coefficient values for the regression 
equations are calculated. The program starts calculating the congelation coefficient value beginning with the first five 
temperature data values of the inverted array, and proceeds towards the last point, adding one point at a time, applying 
the linear least squares fit correlation coefficient equation below: 
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^ (E{T-TJ2 X E(t-t,)^) 

Where: 

= TTie square of the correlation coefficient. 

10 T = The reheat temperature data value. 

= The mean value of the reheat temperature data set. 

t = The corresponding time value. 

= The mean value of the time values. 

15 The calculation of the squares of the correlation coefficients is performed by the subroutine computerrunningarrays 
in the CURVE3 software. The coefficients of each of the linear regression equations are also calculated and are stored 
in the an^ays called aval(1-size%) and bval(1-size%) with the squared congelation coefficients held in the array called 
regg(1 -size%). Because the first four squared correlation coefficients are not computed the first four values of the an^ay 
regg(1-size%) are set to zero. Typically these values should be near unity but it is justifiable to assume that the first lin- 

20 ear regression line should have more than four data points to be significant. This formula is applicable to normally dis- 
tributed data. Since the data under process are the natural logarithms of the initial data it is reasonable to assume that 
the distributions they follow is dose to normal. The program uses these congelation values only for conparison between 
different linear regression equations to determine the most appropriate lit". The algorithm is essentially trying to deter- 
mine the best fit line by emulating the visual method which has been used previously with some success. 

25 The variable rmax is initially set to the value of the squares regression coefficient array value regg(ty%) where 
ty%=5. 

In order to find the maximum squared correlation coefficient the program then differentiates the array regg(1 -size%) 
from the fifth regression coefficient value to the (size%-8)th regression coefficient value using a 3 point backwards dif- 
ference technique, with the differentiated values being stored in the an'ay called diffregg(ty%-(size%-1)). 
30 The algorithm then searches for the last maximum squared concelation coefficient value by testing for the condition 
where the last positive to negative transition occurs in the array diffregg(ty%-(size%-1)). If the above condition occurs 
the array pointer ty1% is set to tiie diffregg(ty%-(size%-1)) array pointer less two (t%-2). The pointer ty1% specifies tiie 
data point number used to determine the value of tiie linear regression coefficients held in the arrays aval(ty1%) and 
bval(ty1%). The variable maxreg is set to the value in the array regg(ty1%). The next test condition is to determine 
35 whetiier tiie gradient of the linear regression line is either positive or negative. If the gradient, i.e. the value of tiie array 
aval(ty1%) is negative, then the curve deconrtposition algorithm will continue to search for the condition where the value 
of the array aval(ty1%) is positive and the value of the array pointer ty1% is set to the new value where this condition 
applies. If the value of array aval(ty1 %} is positive tiien the algoritiim has successfully found the most appropriate linear 
regression line, i.e. the line which fits the natural logarithmic temperature data with the maximum correlation as indi- 
go cated by the squared correlation coefficient value. 

If there are no ti'ansitlons from positive to negative values within the array diffregg(ty%-(size%-8), i.e. there are no 
maxima values of the array regg%(ty%-(size%-1 )) then the process of finding the largest squared regression coefficient 
is implemented by resetting the array pointer t% to ty% artd tiien, by using a FOFVNEXT loop, the curve decomposition 
algorithm searches for the maximum squared correlation coefficient. Again if the value of aray aval(ty1%) is negative 
45 the curve decomposition algorithm will search for the condition where the array ava](t%) is greater than zero and the 
value of the maximum squared correlation coefficient exists. The value of the pointer ty1% is set to the array pointer t% 
and the best regression line has tiie coefficients defined by aval(ty1%) and bval(ty1%). 

After the algorithm has found the most appropriate linear regression equation coefficients which best W the first 
line evaluation of the natural logarithmic temperature data points, the algorithm then subtracts tiie exponential conpo- 
50 nent function exp(A).(-at) from the existing logarithmic temperature data. This subtraction is performed by the subrou- 
tine subtractparameters. The temporary variables tavl and tiovl are set to the values of tiie coefficients of the linear 
regression equation and tiien the subroutine subtractparameters is called. The "new" logarithmic data variable newlg is 
found by carrying out the following. 

55 

newlg =» lgrhv(n%) - (EXP (tbvl)* -tavl* timeval% (n%) ) 
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The lgrhv(n%) array is then set to value of newig, the temporary values tavl and tbvl are reset to zero and the first 
value of the array lgrhv(l-size-%) is set to zero ready for the next curve decomposition. The first curve strip or decom- 
position is now conpleted with the parameter A being the value of the array bval(ty1%) and the a parameter being the 
value of the array aval(ty1%). 

5 The second regression fit is achieved using a similar process working on the new logarithmic curve data with the 
data range set from the data point number (ty1%+1) to (size%-1). Again the squared correlation coefficient monitoring 
is used by recalculating the new values for the arrays regg((ty1%+1)-(size%-1)), aval((ty1%+1)-(slze%-1) and 
bval((ty1%+1)-(size%-1). By finding the mosX appropriate linear regression fit the array pointer ty2% is obtained and the 
coefficients of the linear regression line are returned and the new parameters B and b are again subtracted from the 

10 previous logarithmic data set to form the next set of logarithmic data. The third regression fit is achieved using the same 
process. 

The same curve stripping algorithm was also implemented in the program PR0JECT1 and the graphical represen- 
tations of the reheat temperature data arrays, logarithmic data arrays and the squared correlation coefficients array can 
be displayed by the user. The computed parameters A,a,B,c,C,c together with the maximum squared correlations 
15 are also displayed. To display the graphs the user can use the option (F7-DISPI-AY REHEAT CURVE GRAPH) in the 
main menu of the PR0JECT1 software. By moving the cross-hairs over the first thermographic image the user can dis- 
play the reheat tenrperature data et of the individual pixel by pressing the RETURN key. The parameters A,a,B)C,C,c 
and C,c were then used to reconstruct the reheat tenrperature curve using the following equation. 

20 

Tr(t) = minval + exp(l) + (exp(A) x 6xp(-at)) 

- (exp(B) X exp(-bt)) - (exp(C) x exp(-ct)) 

25 (Note the exponential component time constants are the same values used in the mathematical model but the scal- 
ing coefficients A,B,C used in the imaging are the logarithms of the scaling coefficients used in the nrathematical nrxxJel 
described earlier.) 

When this temperature function is plotted on the display against the experimental data the differences between the 
experimental data and the calculated model data can be graphically seen by the user. The mean value of the differ- 
30 ences (mean value of en-or) was calculated using the fbltowing. (The same time values were used in calculating the 
"predicted" tenperature values.) 

mean_value_of_error« (E (y-x) ) 

35 N 



where, x = the "predicted" data, y = the experimental data and N = the number of data points. 

40 

The mean value of error was displayed to give an indication of the "fit" between the experimental data and the pre- 
dicted data. When the predicted results were calculated and displayed the results were sometimes disappointing and 
gave a poor visual 'lit" whilst other reheat temperature curves gave excellent results. The mean value of en-or results 
were actually a poor indicator of curve reconstruction lit", giving very low values with a poor visual "fit" and in many 

45 cases giving high values of error for curves which fitted very well indeed if the "offset component" is removed. This hap- 
pens when the deviation of the calculated "predicted" curve evolves in a parallel manner to the experimental data. Two 
perfectly parallel similar curves have a correlation coefficient of 1 , but a poor absolute mean value of error (differences). 
Other errors also became apparent. The "predicted" curve when reconstructed showed an initial excellent lit" at the 
beginning of the experimental curve only to deviate badly as the "predicted" curve progresses to the end time values. 

50 Because exponential functions rapidly change as they progress any error in the coefficients of the exponential compo- 
nent functions will quickly become apparent. The errors in the coefficients of the exponential component functions are 
therefore important and can easily arise if for ir^tance a single temperature data point with a considerable error gives 
rise to a wrong set of exponential coefficients during tiie first curve deconrposition when there are only a few tenpera- 
ture values being used to calculate the squared correlation coefficients. Further work may be needed to produce an 

55 adaptive technique to "lock in" the two curves using a feedback technique to obtain even better accuracy in determining 
the parametric values. 
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Parametric File Image Display 

The final stage of the software is a program which displays the parameters produced during the curve decomposi- 
tion. As noted before the six parameters A,a,B,c and C,c are the coefficients of the mathematic model: 

5 

Tr(t)=iiiiiival + Aexp"**+Bexp-**^+Cexp-*="^ 



10 

with minval being the first temperature value of the reheat temperature curve. The mathematical model above was cho- 
sen because of the underlying theoretical model of the thermodynamics of blood perfusion within the dermis. The equa- 
tion has two variables, that of temperature and time and represents the temperature variation of an individual pixel with 
respect to time. The aim of the software is to produce a visual representation within a spatial context of the regional var- 

15 iations of the individual parameters. 

The program RESULTS1 performs this task. The RESULTS1 program is an easy to use program using a simple 
menu system allowing the user to display the computer parameters from the curve decomposition. The program 
RESULTS1 can be run by using the appropriate option (F11-RUN THE PARAMETRIC If^AGE POST PROCESSOR 
DISPLAY PROGRAIW) within the main menu of the PR0JECT1 program. 

20 The user is prompted to enter the path and filename of the FILENAME.DAT parametric image file which was cre- 
ated during the curve decomposition of the reheat curve temperature data. The program RESULTSI is referred to as a 
graphics post processor. A routine for the detection of outliers has been inrplemented. In order to scale the parameters 
con^ectly and to assign the appropriate colours to the differerrt discrimination levels, a range of values for each param- 
eter has to be determined and a statistical approach was taken. The program opens the FILENAME.DATf ile then reads 

25 the text header infornnation which is then displayed to the user. The program then begins to read each parameter in the 
parametric data file and calculates the mean value and standard deviation of each parameter. By using a multiple of the 
standard deviation plus or minus the mean value, a range of values for each parameter can be set. In this case plus or 
minus three standard deviations were used. Any parameter value outside this range is then considered to be an outlier 
and is discarded. The parameter value range together with the parameter mean value and standard deviation are then 

30 displayed to the user. 

Therefore the following was implemented: 

The mean value and standard deviation of each parameter were computed. 

The range between the mean value of each parameter plus and minus three standard deviations was defined. 
35 Every parameter value that was outskle the connputed range was considered to be an outlier. 

The user can then display the parameters using the easy to use menu system which was implemented within the 
parametric imaging software. In order to visualise the regional parameter variations a three dimensional display algo- 
rithm has also been Implemented allowing the user to display graphically the parameters in an isometric form. By press- 
40 ing F1 0 the user can display the parameters in an isometric form and adjustments to the scaling and display angle can 
be easily changed. Also included is a low-pass filter algorithm to smooth out some of the "noise" from the images gen- 
erated. 

The images produced by the software can be seen in Figs. 13 to 15 (showing normal skin data) and Rgs. 17 to 19 
(showing tuberculin reaction skin data). 

45 

Dynamic Thermographic image Processing Documentation File 

The following is a summary of the files used in the thermographic image processing software method according to 
the invention and referred to in the description of the method, with a short description of each file. 

so 

CONVERT1 The AGEMA conversion utility source code which aeates ASCII output data files written in TufboC++. 
CON VERT2 The AGEMA conversion utility source code which aeates binary output data files written in TurtwC-H-. 
PR0JECT1 The thermographic image processing program written in QBASIC. 
RESULTSI The parametric image post processor program written in QBASIC. 
55 CURVE3 The curve decomposition source program written in Turtx)BASIC. 

Thus the invention provides a method of dynamic imaging of the blood perfusion rate in skin which, because of the 
speed at which the images can be created, is useful as a diagnostic tool for the medical profession. 
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The methcxi is able to illustrate the changes in skin temperature rather than the skin temperature alone. These 
changes are directly related to blood perfusion rates in the deeper layers of the dermis having diameter >50M.m. 

The veterinary applications for Infrared thermography are numerous also. Not only can it help veterinary surgeons 
make objective clinical diagnoses, but It can be used to nrtonitor the progress of n^ny different physiological disorders 
5 and gauge their response to various courses of treatment. It can immediately help with the identification of a whole 
range of problems, including: muscle disorders and secondary site of disease In lameness problems, conditions before 
synrtptoms show, eg tendons and joint disease, whether spllrrts are active or inactive, back/spinal problems and Identi- 
fying the difference between recent and long standing disorders, nerve injuries and circulatory disorders. 

The curve stripping of multi- exponential curve decomposition is common to clinical studies. The underlying Idea 
70 behind this modelling is the representation of a system with a finite number of component parts, called conpartments 
or components. Each corrpartment contributes to the system according to its weighting factor, but has different origin, 
quality, and characteristics. 

The temperature of the forearm's skin is the system under measurement. Previous studies indicated that skin tem- 
perature is a two-compartmental system. The method of the invention assumes that a third component could also fit the 
15 reheat curves. 

Ultimately, 38.5% of the tuberculin reaction skin data (TRSD) and 55.9% of the normal skin data (NSD) revealed 
three cooponents. Nevertheless, the two first components exist beyond any doubt, and the existence of a third corrpo- 
nent, Is indicated for all the data. The physiological nature of the two components can be described as follows: 

20 • a first component largely dependent on one or two conpartments of blood flow 

a second component largely dependent on non-perfusion heat exchange mechanisms, including conduction 

Figs. 13 to 15 and 17 to 19 Illustrate the constructed coefficient images for the normal and the tuberculin reaction 
skin data respectively One control thermographic Image (acquired before the application of the dynamic thernrwgraphtc 

25 procedure) is also displayed as Figs. 12 and 16 respectively at the start of each set of Images. These two control 
Images are both normalized from zero to one. The normalisation Is done in each image's individual range (minimum- 
maximum values) and for this reason conparison between them In absolute temperature terrrs Is not applicable. 

On the other hand, comparisons between the two coefficient sets are more meaningful. The control images are dis- 
played In order to be compared with the corresponding coefficient images. 

30 For each corrponent the "capital letter" image (A, B, C) Is its weighting factor, and the other image (snnall letter a, 
b. c) can be interpreted as the one related to perfusion rates. The Identification of which conponent corresponds to heat 
exchange due to blood perfusion and which to heat exchange due to other mechanisms is difficult at this stage on the 
basis of only two images. 

After long consideration and observation of the images it was the first component that seems to be the corrponent 
35 (largely) dependent on blood circulation (A, a). This approach Illustrates the principle of higher flow (a) in small com- 
partments (A) and vice versa. 

An evidence that this evaluation Is correct is the clear imaging of a vein that appeared in some parametric images 
and In the corresponding control Images. The "A" value of the vein was low ("cold'l and the a value was high ("hot"), 
leading to the physiological interpretation of high flow in small conpartment. 
40 The normal skin data Images did not show big differences compared to the control inrage. High flow In small com- 
partments, and a heterogeneity level was obvious. A possible boundary effect was also displayed in A. a and b images, 
but not in the B image. 

The tuberculin reaction skin data Images were much more interesting. They illustrated clearly a high level of heter- 
ogeneity in red blood cell distritxjtion (A) and velocities (a). There was considerably higher resolution of regional 
45 changes in blood circulation compared with the control image. 

These images Indicate the existence of very fine regulatory mechanisms controlling the blood flow in the deeper 
layers of the skin. The principle of high flow in snnall conpartments is clear. 

Dynamic thermography of the tuberculin reaction, followed by multi-exponential curve decomposition, indicates 
blood perfusion differences in a regional scale. Much better contrast In blood flow values is obtained compared to the 
50 static thermographic measurements. 

The images illustrate the size of the conpartments of blood flow together with the rate of flow through them. 

Although particular software solutions are desalbed above, it is to be understood that the invention is not to be lim- 
ited to these particular software solutions, but is instead defined by the scope of the claims. 

Similarly the apparatus of the invention should not be limited to those apparatus conponents specifically described 
55 above, but should include any infra-red detecting means which may be used to obtain periodically a set of thermo- 
graphic data corresponding to skin temperature on an area of skin, including bolometers, thermoelectric detectors, 
scanners and pyroelectric vidicons; any data storage means which may be coupled to the detecting means and can 
store the thermographic data; any computational means which can be programmed to calculate the coefficients of an 
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exponential fit equation, including separate PCs and dedicated microprocessor control units; and any image display 
means adapted to display the coefficients of the exponential fit equation in graphical form, including printers, monitors 
and TV screens. 

These and other modifications and improvements can be incorporated without departing from the scope of the 
invention. 

Claims 

1 . A method of measuring blood perfusion in human or animal skin, comprising the steps of: 



(a) applying a cooling means to an area of skin; 

(b) removing the cooling means from the area of skin; 

(c) using an infra-red detector to obtain a first set of thermographic data TDi corresponding to skin temperature 
at a plurality of points on the area of skin at time t^ ; 

15 (d) repeating step (c) at least once to obtain a total of n sets of thermographic data TDi to TD^^ corresponding 

to times t^ to tn, where n is at least 2; 

(e) calculating the coefficients of an exponential fit equation for the variation with time of the thermographic 
data TD^ to TD^ at each of the plurality of points; and 

(f) displaying the coefficients of the exponential fit equation in graphical form. 

20 

2. A method according to Claim 1 , wherein the coefficients of the exponential fit equation are calculated using a curve 
stripping technique. 

3. A method according to Claims 1 or 2, wherein n is at least 4, preferably at least 10. 

4. A method according to any preceding Claim, wherein the coefficients of the exponential fit equation are displayed 
using colour mapping. 

5. A method according to any preceding Claim, wherein the cooling means is a cold object which is placed on the skin. 

6. A method according to any preceding Claim, wherein the method further includes the step of placing a reference 
object on the area of skin so that the thermographic data TD^ is obtained at the same points on the area of skin at 
each successive repeat of step (c). 

35 7. A method according to Claim 6. wherein the reference object is a reflective strip. 

8. A method of processing and displaying data relating to blood perfusion in human or animal skin, comprising the 
steps of : 

40 (a) generating a series of n thermographic image matrices IM^ to IM^ corresponding to times ti to tn, each ther- 

mographic image matrix comprising temperature data relating to a plurality of points on an area of skin; 

(b) creating a filter matrix corresponding to a region of interest selected from the area of skin; 

(c) applying the filter matrix to each thermographic image matrix IM^ to IMr, to obtain an cropped thermographic 
image matrix IMCR^ to IMCRn having the value zero for points outside the region of interest; 

45 (d) constructing an array AR^ for each point m in the region of interest, the array containing temperature data 

Tfn at that point corresponding to times t^ to tp; 

(e) for each array AR^ performing a regression analysis to calculate the coefficients of the exponential compo- 
nents of the best fit equation 

50 T^=T^,, + Ae«^Be"+Ce^ 

where A, a, B, b, C, c are coefficients and Tf^jn is the minimum temperature data in each array; 

(f) constructing a coefficient matrix (Z^, Za etc.) having as elements a parameter derived from the coefficients 
55 (A, a etc.) calculated in step (e) for each point m in the region of interest; and 

(g) displaying the coefficient matrix as an image. 

9. A method according to Claim 8, wherein the parameter in step (f) is equal to one of the following: 
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(i) one of the coefficients (A. a etc.) calculated in step (e); 

(il) the sum of the coefficients A, B and C calculated in step (e). representing the total tissue volume: 
(iii)the sum of the coefficient products (aA + bB + cC) divided by the sum of the coefficients (A + B C), repre- 
senting the total flow. 

5 

10. A method according to Claim 8 or 9, wherein the region of interest corresponds to an area of skin which has previ- 
ously been cooled. 

1 1 . A method according to any one of Claims 8 to 1 0, wherein the coefficients of the exponential components are Cai- 
ro culated using a curve stripping technique. 

1 2. A method according to any one of Claims 8 to 11 , wherein the minimum temperature data T^in in each array is sub- 
tracted from each temperature data in the array before calculation of the coefficients. 

15 13. A method according to any one of Claims 8 to 12, wherein the data corresponding to a first interval of time after 
cooling is discarded before calculating coefficients. 

1 4. A method according to Claim 1 3, wherein the first i nterval of time is preferably less than 1 0 seconds, more prefer- 
ably less than 5 seconds, most preferably between 3 and 4 seconds. 



20 



15. An apparatus for measuring blood perfusion in human or animal skin, comprising: 



(a) an infra-red detecting means adapted to obtain periodically a set of thermographic data (TD^ to TDn) cor- 
responding to skin temperature at a plurality of points on an area of skin; 
25 (b) data storage means adapted to store each set of thermographic data (TD^ to TDj^\ 

(c) computational means for calculating the coefficients of an e^qponential fit equation for the variation with time 
of the thermographic data (TD^ to TDn) at each of the plurality of points; and 

(d) image display means adapted to display the coefficients of the exponential fit equation in graphical form. 

30 



35 



40 



45 



50 



55 



17 



EP0 885 587A1 




Fig. 1 




Fig. 2 



18 



EP0885587A1 




0 ' ' 

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 

WAVELENGTH {fim) 



Fig. 4 



19 



EP 0 885 587 A1 




EP 0 885 587 A1 




Fig. 8 
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